1 Traitement de données spatiales

1.1 Préparer son environement de travail

Créez un projet RStudio.
Dans RStudio créez un projet et nommé le exercice. Il s’agit d’une bonne pratique qui vous facilitera la tâche. Cela améliore l’organisation et la portabilité de votre travail.
Dans le projet, créez un dossier data.
Téléchargez les données :
- Un fichier contenant les polygones des communes d’Occitanie au format gpkg est disponible ici : Occitanie.gpkg. Il faut télécharger le fichier puis le placer dans le dossier data que vous venez de créer.
- Un fichier contenant des données socio-démographiques sur les communes de France est disponible ici : base_cc_comparateur.xls. Il faut le télécharger et le mettre également dans le même dossier data.
Créez un fichier exo1.R et enregistrez le à la racine du projet.

Votre projet RStudio doit être structuré de la manière suivante :

├── data
│   └── base_cc_comparateur.xls
│   └── Occitanie.gpkg
├── exo1.R
└── exercice.Rproj

Vous devez enregistrer votre script régulièrement dans le fichier exo1.R.

Source du fichier Occitanie.gpkg : ADMIN EXPRESS COG édition 2019, IGN
Source du fichier base_cc_comparateur.xls : Base comparateur de territoires 2019, INSEE

1.2 Charger une couche géographique

Chargez la couche géographique des communes d’Occitanie en utilisant le package sf et la fonction st_read().
Vérifiez le système de projection avec st_crs(). S’agit-il de données projetée ? Si ce n’est pas le cas, transformez la couche des communes dans la projection française (Lambert 93, EPSG : 2154) avec la fonction st_transform().

library(sf)

## Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0

occ_raw <- st_read(dsn = "data/Occitanie.gpkg", stringsAsFactors = FALSE)

## Reading layer `Occitanie' from data source `/home/tim/Documents/prz/carto_avec_r_exo/data/Occitanie.gpkg' using driver `GPKG'
## Simple feature collection with 4454 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -0.3271723 ymin: 42.33292 xmax: 4.845565 ymax: 45.04557
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

st_crs(occ_raw)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

occ <- st_transform(x = occ_raw, crs = 2154)

1.3 Réaliser une sélection par attributs

Séléctionnez toutes les communes de la Haute-Garonne.
Enregistrez votre séléction dans un nouvel objet nommé com31
Afficher la nouvelle couche (en utilisant notament st_geometry())

Pour connaitre la liste de tous les noms ou code de région on peut utiliser la fonction unique().

unique(occ$INSEE_DEP)

##  [1] "32" "31" "82" "12" "81" "30" "11" "65" "46" "09" "34" "66" "48"

com31 <- occ[occ$INSEE_DEP == "31", ]
plot(st_geometry(com31))

1.4 Fusionner des entités

Fusionnez les communes de la région en un seul polygone (reg76), utilisez la fonction st_union ().
Créez la couche géograpique des départements de la région (dep76), utilisez la fonction aggregate() pour regrouper les polygones et calculer les sommes des populations communales.
Affichez les résultats (communes, départements et région)

reg76 <- st_union(occ)
dep76 <- aggregate(occ["POPULATION"], by = list(occ$INSEE_DEP), sum)
plot(st_geometry(occ), col = "lightblue1", border = "white", lwd = .5)
plot(st_geometry(dep76), col = NA, border = "lightblue2", lwd = 1, add = TRUE)
plot(reg76, col = NA, border = "lightblue3", lwd = 2, add = TRUE)

1.5 Créer une zone tampon

Créez une zone tampon d’une distance de 10 km autour des limites de la commune de Toulouse avec st_buffer(). Quel est le code INSEE de la commune de Toulouse?

toulouse <- com31[com31$INSEE_COM == "31555",]
toulouse_buffer <- st_buffer(toulouse, 10000)
plot(st_geometry(com31), lwd = .5)
plot(st_geometry(toulouse), col = "darkblue", add = TRUE)
plot(st_geometry(toulouse_buffer), col = NA, lwd = 3, border = "red", 
     add = TRUE)

1.6 Réaliser un sélection par localisation

Déterminez quelles communes de la Haute-Garonne intersectent le buffer créé.

Utilisez la fonction st_intersects().
Inserez directement le resultat dans une variable nommée in_buffer dans l’objet com31.

com31$in_buffer <- st_intersects(x = com31, toulouse_buffer, sparse = FALSE)
head(com31, 2)

## Simple feature collection with 2 features and 12 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 504779 ymin: 6198834 xmax: 565151 ymax: 6270449
## epsg (SRID):    2154
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
##                         ID         STATUT INSEE_COM INSEE_ARR INSEE_DEP
## 4 COMMUNE_0000000009762857 Commune simple     31046         2        31
## 5 COMMUNE_0000000009760450 Commune simple     31547         1        31
##   INSEE_REG CODE_EPCI NOM_COM_M POPULATION TYPE NOM_COM
## 4        76 200072635     BAREN         12  COM   Baren
## 5        76 200068641   SEYSSES       8787  COM Seysses
##                             geom in_buffer
## 4 MULTIPOLYGON (((504930 6199...     FALSE
## 5 MULTIPOLYGON (((556904 6267...      TRUE

1.7 Afficher des couches géographiques

Affichez/superposez toutes les couches géographiques créees :

Les communes de la Haute-Garonne
Les communes intersectées par le buffer
Toulouse
Le buffer autour de toulouse
Les départements d’Occitanie
Les limites de l’Occitanie

Jouez sur les styles pour les différencier…

plot(reg76, col = NA, border = "grey50", lwd = 2, bg = "lightyellow")
plot(st_geometry(com31), col = "#aec8f2", border = "white", lwd = .5, 
     add = TRUE)
plot(st_geometry(com31[com31$in_buffer == TRUE, ]),col = "red", lwd = .5, 
     border = "white", add = TRUE)
plot(st_geometry(toulouse), col = "darkblue", border = "black", lwd = 1,
     add = TRUE)
plot(st_geometry(toulouse_buffer), col = NA, border = "black", lwd = 2, 
     lty = 2, add = TRUE)
plot(st_geometry(dep76), col = NA, border = "grey50", add = TRUE)

1.8 Créer une couche de points

Créer un objet sf (points) contenant la localisation de la préfécture de région à Toulouse.

Vous pouvez récupérer la longitude et la latitude sur Google Map ou OSM.
Créer un point avec st_points(), puis un objet sfc avec avec st_sfc() et finalement un objet sf avec st_sf().

# Sur OSM : https://www.openstreetmap.org/search?whereami=1&query=43.59825%2C1.45047#map=19/43.59825/1.45047
# Position de la préfecture 43.59825, 1.45047
library(sf)
pref_pt <- st_point(c(1.45047, 43.59825))
pref_sfc <- st_sfc(pref_pt, crs = (4326))
pref <- st_sf(name = "Préfecture", geometry = pref_sfc)
pref

## Simple feature collection with 1 feature and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 1.45047 ymin: 43.59825 xmax: 1.45047 ymax: 43.59825
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##         name                 geometry
## 1 Préfecture POINT (1.45047 43.59825)

1.9 Calculer des distances

Calculez une matrice de distance entre la préfecture et les centroïdes des communes du département Haute-Garonne.
Ajouter cette distance dans une nouvelle colonne dist_pref dans com31. Pour cela, utiliser les fonction st_centroid() et st_distance().
N’oubliez pas de vérifier les projections utilisées avec st_crs(), au besoin modifiez les avec st_transform().

st_crs(pref)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

st_crs(com31)

## Coordinate Reference System:
##   EPSG: 2154 
##   proj4string: "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

identical(st_crs(pref), st_crs(com31))

## [1] FALSE

pref <- st_transform(pref, st_crs(com31))
identical(st_crs(pref), st_crs(com31))

## [1] TRUE

# Centroides des communes
com31_centro <- st_centroid(st_geometry(com31))
com31$dist_pref <- st_distance(com31_centro, pref)
head(com31, 2)

## Simple feature collection with 2 features and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 504779 ymin: 6198834 xmax: 565151 ymax: 6270449
## epsg (SRID):    2154
## proj4string:    +proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
##                         ID         STATUT INSEE_COM INSEE_ARR INSEE_DEP
## 4 COMMUNE_0000000009762857 Commune simple     31046         2        31
## 5 COMMUNE_0000000009760450 Commune simple     31547         1        31
##   INSEE_REG CODE_EPCI NOM_COM_M POPULATION TYPE NOM_COM
## 4        76 200072635     BAREN         12  COM   Baren
## 5        76 200068641   SEYSSES       8787  COM Seysses
##                             geom in_buffer     dist_pref
## 4 MULTIPOLYGON (((504930 6199...     FALSE 104829.70 [m]
## 5 MULTIPOLYGON (((556904 6267...      TRUE  17211.48 [m]

1.10 Cartographier une variable

Cartographiez distance de chaque commune du département à la préfecture avec la fonction plot().

# Le plus simple
plot(com31["dist_pref"])

# Avec un peu de paramétrage
plot(com31["dist_pref"], main = "Distance à la préfecture (en mètres)", 
     pal = hcl.colors(13,"Turku", rev = TRUE), border = NA,
     key.pos = 1,key.width = .15, key.length = .75, graticule = TRUE,
     reset = FALSE)
plot(st_geometry(pref), pch = 20, col = "red", add = TRUE )

1.11 Extraire des données OpenStreetMap

Sélectionnez les communes d’un autre département dans un nouvel objet nommé com+codeDuDépartement* (com46 par exemple).
Avec le package osmdata extrayez l’ensemble des restaurants présents dans le département.
Affichez les communes qui ne contiennent pas de restaurants.
Pour ces communes calculez la distance au restaurant le plus proche.

com46 <- occ[occ$INSEE_DEP == "46", ]
plot(st_geometry(com46))

library(osmdata)

## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

# Définition d'une bounding box
q <- opq(bbox=st_bbox(st_transform(com46,4326)))
# Extraction des restaurants
res <- add_osm_feature(opq = q, key = 'amenity', value = "restaurant")
res.sf <- osmdata_sf(res)
res.sf.pts  <- res.sf$osm_points[!is.na(res.sf$osm_points$amenity),]
resto <- st_transform(res.sf.pts, st_crs(com46))
resto <- st_intersection(resto, st_geometry(com46))

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

# Affichage des restaurants
plot(st_geometry(com46), col="darkseagreen3", border="darkseagreen4")
plot(st_geometry(resto), add=TRUE, pch=20, col = "#330A5FFF", cex = 0.5)

# Compter les restaurants par communes
inter <- st_intersects(x = com46, y = resto)
com46$nresto <- sapply(inter, length)

# communes sans restos
com46noresto <- com46[com46$nresto==0, ]
plot(st_geometry(com46), col="darkseagreen3", border="darkseagreen4")
plot(st_geometry(com46noresto), col = 'red', add = TRUE)
plot(st_geometry(resto), add=TRUE, pch=20, col = "blue", cex = 1)

# index du restaurant le plus proche
index <- st_nearest_feature(x = st_centroid(com46noresto), 
                            y = resto)

## Warning in st_centroid.sf(com46noresto): st_centroid assumes attributes are
## constant over geometries of x

# distance au plus proche
com46noresto$dresto <- st_distance(x = st_centroid(com46noresto), 
                                   y = resto[index, ], 
                                   by_element = TRUE)

## Warning in st_centroid.sf(com46noresto): st_centroid assumes attributes are
## constant over geometries of x

# Affichage de la carte
plot(com46noresto['dresto'], reset = F, 
     main = "Distance au restaurant le plus proche")
plot(st_geometry(com46), col=NA, add= TRUE)
plot(st_geometry(resto), add=TRUE, pch=20, col = "red", cex = 1)

2 Cartographie thématique

2.1 Préparer son environement de travail

Dans le même projet RStudio créez un deuxième script nommé exo2.R.
Votre projet RStudio doit maintenant être structuré de la manière suivante :

├── data
│   └── base_cc_comparateur.xls
│   └── Occitanie.gpkg
├── exo1.R
├── exo2.R
└── exercice.Rproj

2.2 Charger des données statistiques

Charger la couche géographique des communes d’Occitanie (data/Occitanie.gpkg). Utilisez le package sf et la fonction st_read() pour importer les données.
Vérifiez le système de projection avec st_crs().
Charger les fichier de données data/base_cc_comparateur.xls fournie par l’INSEE. Utilisez le package readxl et la fonction read_excel() pour ouvrir la table de données correctement. Importer la table dans un objet nommé occ_df.

2.3 Réaliser une sélection par attributs

Créez un nouvel objet sf à partir d’une séléction par attribut`:

Séléctionnez toutes les communes d’un seul département (ne choisissez pas la Haute-Garonne).
Enregistrez votre séléction dans un nouvel objet nommé com+codeDuDépartement* (com46 par exemple)
Afficher la nouvelle couche (plot(st_geometry(...)).

Pour connaitre la liste de tous les noms ou code de région voir la fonction unique()

2.4 Réaliser une jointure

Joignez la table de données INSEE avec la couche gégographique des communes de votre département séléctionné. Utilisez la fonction merge(). Quel est l’identifiant commun entre la table géo et la table attributaire?

2.5 Réaliser des cartes thématiques

Réalisez huit cartes thématiques finalisées.

Chacune des cartes devra comprendre les éléments d’habillage et de mise en page nécessaires leur compréhension (titre, légende, sources…)

Parmi ces huits cartes deux échelles différentes (échelle de la région, échelle d’un département) et deux découpages territoriaux différents (les communes, les départements) doivent être utilisés.

Vous devez réaliser les huit types de cartes suivants :

Une carte de localisation (avec carton de localisation si possible, des labels, legende…);
Une carte représentant une variable quantitative absolue (un stock);
Une carte représentant une variable quantitative relative (un ratio);
Une carte combiant une variable quantitative absolue et une variable quantitative relative (un stock et un ratio);
Une carte représentant une variable qualitive;
Une carte utilisant une anamorphose. Il est possible de combiner les anamorphoses avec d’autres variables (qualitative ou quantitative relative);
Une carte interactive;
Une carte utilisant une grille régulière ou représentant des discontinuités.

Exemple :

Une carte en symboles proportionels représentant la population des départements de la région (échelle régionale, découpage départementale, carte de stock.).
Une carte choroplète représentant l’évolution de la population des communes d’un département (échelle départementale, découpage communale, carte de ratio).

Nous vous conseillons d’utiliser le package cartography pour l’ensemble de vos réalisations. Néanmoins l’utilisation d’autres packages comme ggplot2 ou tmap est possible.

2.6 Évaluation

Vous devrez nous rendre pour le vendredi 7 février 2020 un script permettant de construire ces huit cartes dans un fichier dont le nom est constitué de votre prénom et de votre nom (Prenom_Nom.R).
Votre script doit pouvoir fonctionner une fois inséré dans un projet RStudio structuré de la manière suivante :

├── data
│   ├── base_cc_comparateur.xls
│   └── Occitanie.gpkg
├── Prenom_Nom.R
└── rendu.Rproj

Vous serez évalués de la manière suivante :

2 points par cartes (choix sémiologiques, traitement des données en amont…),
4 points pour la qualité et la présentation du code (reproductibilité, commentaires, clareté).

UMS RIATE, CNRS↩
FR CIST, CNRS↩

Cartographie avec R

Timothée Giraud¹ & Hugues Pecout²

14 janvier 2020

1 Traitement de données spatiales

1.1 Préparer son environement de travail

1.2 Charger une couche géographique

1.3 Réaliser une sélection par attributs

1.4 Fusionner des entités

1.5 Créer une zone tampon

1.6 Réaliser un sélection par localisation

1.7 Afficher des couches géographiques

1.8 Créer une couche de points

1.9 Calculer des distances

1.10 Cartographier une variable

1.11 Extraire des données OpenStreetMap

2 Cartographie thématique

2.1 Préparer son environement de travail

2.2 Charger des données statistiques

2.3 Réaliser une sélection par attributs

2.4 Réaliser une jointure

2.5 Réaliser des cartes thématiques

2.6 Évaluation

Cartographie avec R

Timothée Giraud1 & Hugues Pecout2

14 janvier 2020

1 Traitement de données spatiales

1.1 Préparer son environement de travail

1.2 Charger une couche géographique

1.3 Réaliser une sélection par attributs

1.4 Fusionner des entités

1.5 Créer une zone tampon

1.6 Réaliser un sélection par localisation

1.7 Afficher des couches géographiques

1.8 Créer une couche de points

1.9 Calculer des distances

1.10 Cartographier une variable

1.11 Extraire des données OpenStreetMap

2 Cartographie thématique

2.1 Préparer son environement de travail

2.2 Charger des données statistiques

2.3 Réaliser une sélection par attributs

2.4 Réaliser une jointure

2.5 Réaliser des cartes thématiques

2.6 Évaluation

Timothée Giraud¹ & Hugues Pecout²